#####################################################################
# Check sensitivity to hidden biases
#####################################################################

# This code replicates the results summarized in Appendix B: Unmeasured Covariates

rm(list=ls())

library(Hmisc)
library(exactRankTests)
library(ggplot2)
library(stargazer)
library(sensitivitymw)
library(sensitivitymv)
library(lmtest)
library(sandwich)
library(Cairo)
library(gridExtra)

# function for sensitivity analysis
sens.analysis.signedrank = function(diff, Gamma){
  rk = rank(abs(diff))
  s1 = 1*(diff>0)
  s2 = 1*(diff<0)
  W = sum(s1*rk)
  Eplus = sum((s1+s2)*rk*Gamma)/(1+Gamma)
  Eminus = sum((s1+s2)*rk)/(1+Gamma)
  V = sum((s1+s2)*rk*rk*Gamma)/((1+Gamma)^2)
  Dplus = (W-Eplus)/sqrt(V)
  Dminus = (W-Eminus)/sqrt(V)
  list(lowerbound = 1-pnorm(Dminus), upperbound = 1-pnorm(Dplus))
}

# load function that does clustered SEs
vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<40){
    warning("Fewer than 40 clusters, variances may be unreliable")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

##############################################
# Read data and prepare basic output
##############################################

# Generate log  
#sink("006_unmeasured_covariates.txt")

# Read data
load("001_match_ego.Rdata")
d_ego = d
d_match_ego = d_match
d = NULL
d_match = NULL

# Read data 
load("002_match_socio.Rdata")
d_socio = d
d_match_socio = d_match
d = NULL
d_match = NULL

########################### 
# Sensitivity analysis
###########################

# Egotropic
# Sensitivity Analysis for Signed Rank Statistic
treatment_ego = d_match_ego$outcome[d_match_ego$t_ind==1]
control_ego = d_match_ego$outcome[d_match_ego$t_ind==0]

diff_ego = control_ego-treatment_ego # inverse since it is a one side test

ego_sd_100 = round(sens.analysis.signedrank(diff_ego, 1.00)$upperbound, digits=5)
ego_sd_105 = round(sens.analysis.signedrank(diff_ego, 1.05)$upperbound, digits=5) 
ego_sd_110 = round(sens.analysis.signedrank(diff_ego, 1.10)$upperbound, digits=5) 
ego_sd_115 = round(sens.analysis.signedrank(diff_ego, 1.15)$upperbound, digits=5) 
ego_sd_120 = round(sens.analysis.signedrank(diff_ego, 1.20)$upperbound, digits=5)
ego_sd_125 = round(sens.analysis.signedrank(diff_ego, 1.25)$upperbound, digits=5)
ego_sd_130 = round(sens.analysis.signedrank(diff_ego, 1.30)$upperbound, digits=5)
ego_sd_131 = round(sens.analysis.signedrank(diff_ego, 1.31)$upperbound, digits=5)
ego_sd_132 = round(sens.analysis.signedrank(diff_ego, 1.32)$upperbound, digits=5)
ego_sd_133 = round(sens.analysis.signedrank(diff_ego, 1.33)$upperbound, digits=5)
ego_sd_134 = round(sens.analysis.signedrank(diff_ego, 1.34)$upperbound, digits=5)
ego_sd_135 = round(sens.analysis.signedrank(diff_ego, 1.35)$upperbound, digits=5)  # sensitive

# Sociotropic
# Sensitivity Analysis for Signed Rank Statistic 
treatment_socio = d_match_socio$outcome[d_match_socio$t_ind==1]
control_socio = d_match_socio$outcome[d_match_socio$t_ind==0]

diff_socio = treatment_socio - control_socio # 

socio_sd_100 = round(sens.analysis.signedrank(diff_socio, 1.00)$upperbound, digits=5) # sensitive
socio_sd_105 = round(sens.analysis.signedrank(diff_socio, 1.05)$upperbound, digits=5)
socio_sd_110 = round(sens.analysis.signedrank(diff_socio, 1.10)$upperbound, digits=5) 
socio_sd_115 = round(sens.analysis.signedrank(diff_socio, 1.15)$upperbound, digits=5) 
socio_sd_120 = round(sens.analysis.signedrank(diff_socio, 1.20)$upperbound, digits=5) 
socio_sd_125 = round(sens.analysis.signedrank(diff_socio, 1.25)$upperbound, digits=5) 
socio_sd_130 = round(sens.analysis.signedrank(diff_socio, 1.30)$upperbound, digits=5)
socio_sd_131 = round(sens.analysis.signedrank(diff_socio, 1.31)$upperbound, digits=5)
socio_sd_132 = round(sens.analysis.signedrank(diff_socio, 1.32)$upperbound, digits=5)
socio_sd_133 = round(sens.analysis.signedrank(diff_socio, 1.33)$upperbound, digits=5)
socio_sd_134 = round(sens.analysis.signedrank(diff_socio, 1.34)$upperbound, digits=5)
socio_sd_135 = round(sens.analysis.signedrank(diff_socio, 1.35)$upperbound, digits=5)

# Generate sensitivity table 
ego_sen = c(ego_sd_100,ego_sd_110,ego_sd_120,ego_sd_130,ego_sd_131,ego_sd_132,ego_sd_133,ego_sd_134,ego_sd_135)
socio_sen = c(socio_sd_100,socio_sd_110,socio_sd_120,socio_sd_130,socio_sd_131,socio_sd_132,socio_sd_133,socio_sd_134,socio_sd_135)
gamma = c(1.00,1.10,1.20,1.30,1.31,1.32,1.33,1.34,1.35)
sen = data.frame(gamma,ego_sen,socio_sen)
colnames(sen) <- c("Gamma", "Egotropic treatment", "Sociotropic treatment")

stargazer(sen, 
          type = "text", 
          colnames = TRUE,
          summary = FALSE, 
          title="Upper bounds on the one-sided P-value testing the null null hypothesis of no treatment effect, using a Wilcoxon signed rank test statistis",
          digits=3, 
          rownames=FALSE, 
          align=TRUE, 
          float = TRUE,  
          float.env = "table", 
          table.placement = "H")

# Amplification test
# amplify(gamma, lambda)
# lambda: odds of exposure to treatment (top)
# delta: odds of greater outcome (bottom)
amplify(1.32, c(2))

#sink()
